/*
 * File:   latticeSP.cpp
 * Author: gian
 *
 * Created on December 10, 2012, 3:41 PM
 */

#include <math.h>
#include <stdio.h>
#include <stdexcept>
#include <ios>
#include "mt64.h"
#include "latticeSP.hpp"



nodeSP::nodeSP() {
  theta = 0.0;
  omega = 0.0;
  n     = 0;
}

nodeSP::nodeSP(const nodeSP & orig) {
  n     = orig.n;
  theta = orig.theta;
  omega = orig.omega;
}

nodeSP::~nodeSP() {}

    void nodeSP::save(FILE *stream)
    {
        fwrite((void*)&n,sizeof(int),1,stream);
        fwrite((void*)&theta,sizeof(double),1,stream);
        fwrite((void*)&omega,sizeof(double),1,stream);
    }
    void nodeSP::load(FILE *stream)
    {
        fread((void*)&n,sizeof(int),1,stream);
        fread((void*)&theta,sizeof(double),1,stream);
        fread((void*)&omega,sizeof(double),1,stream);        
    }


latticeSP::latticeSP(int nx,
                     int ny,
                     int nz,
                     int n) {
  mNx = nx;
  mNy = ny;
  mNz = nz;
  mN  = n;

  int lattice_size = mNx * mNy * mNz;

  mNodes = new nodeSP[lattice_size];

  set_periodic_boundary();
}

latticeSP::latticeSP(const latticeSP & orig) {
  mNx = orig.mNx;
  mNy = orig.mNy;
  mNz = orig.mNz;
  mN  = orig.mN;

  int lattice_size = mNx * mNy * mNz;

  mNodes = new nodeSP[lattice_size];

  set_periodic_boundary();
}

latticeSP::~latticeSP() {
  delete[] mNodes;
}

void latticeSP::init_random(double sigmaPI) {
  int      lattice_size = mNx * mNy * mNz;
  nodeSP * p            = mNodes;

  for (int i = 0; i < lattice_size; i++) {
    p -> n     = (int) (genrand64_real2() * mN);
    p -> theta = M_PI * genrand64_real2();
    p -> omega = sigmaPI*genrand_normal();
    p++;
  }
  
  
}

void latticeSP::set_periodic_boundary() {}

int latticeSP::save_vtk(const char * filename) {
  FILE * output;
  char   buff[128];
  int    i, k;
  int    n, lattice_size;
  float  val;

  snprintf(buff, 128, "%s.vtk", filename);

  output = fopen(buff, "wb");

  nodeSP * p = mNodes;

  fprintf(output, "# vtk DataFile Version 3.0\n");
  fprintf(output, "Phase values\n");
  fprintf(output, "ASCII\n");
  fprintf(output, "DATASET STRUCTURED_POINTS\n");
  fprintf(output, "DIMENSIONS %d %d %d\n", mNx, mNy, mNz);
  fprintf(output, "ASPECT_RATIO 1 1 1\n");
  fprintf(output, "ORIGIN 0 0 0\n");
  fprintf(output, "POINT_DATA %d\n", mNx * mNy * mNz);

  lattice_size = mNx * mNy * mNz;

  fprintf(output, "SCALARS phase float 1\n");
  fprintf(output, "LOOKUP_TABLE local\n");

  for (i = 0; i < lattice_size; i++) {
    val = p[i].theta / M_PI + p[i].n;

    fprintf(output, "%f\n", val);
  }

  fprintf(output, "LOOKUP_TABLE local 300\n");

  for (i = 1; i <= 100; i++) {
    fprintf(output, "%f 0.0 0.0 1.0", (float) i / 100.0);
  }

  for (i = 1; i <= 100; i++) {
    fprintf(output, "0.0 %f 0.0 1.0", (float) i / 100.0);
  }

  for (i = 1; i <= 100; i++) {
    fprintf(output, "0.0 0.0 %f 1.0", (float) i / 100.0);
  }

  fclose(output);
}
#define DELTA 0.03

struct line {
  line(float x_s,
       float y_s,
       float z_s,
       float x_e,
       float y_e,
       float z_e,
       int   c1,
       int   c2) {
    xs     = x_s;
    ys     = y_s;
    zs     = z_s;
    xe     = x_e;
    ye     = y_e;
    ze     = z_e;
    color1 = c1;
    color2 = c2;
  }

  ;
  int   color1;
  int   color2;
  float xs, ys, zs;
  float xe, ye, ze;
};


std::list<line> line_list;

void latticeSP::measure_link(nodeSP *              n1,
                             nodeSP *              n2,
                             std::vector<double> & colors) {
  double deltafi;
  double fi1 = n1 -> theta;
  double fi2 = n2 -> theta;

  if (n1 -> n == n2 -> n) {
    colors[n1 -> n] += fi2 - fi1;
  } else {
    deltafi = 2.0 * M_PI - fi2 - fi1;

    if (deltafi < M_PI) {
      colors[n1 -> n] += M_PI - fi1;
      colors[n2 -> n] -= M_PI - fi2;
    } else {
      colors[n1 -> n] -= fi1;
      colors[n2 -> n] += fi2;
    }
  }
}

void latticeSP::save(FILE * output)
{
    latticeSP & L = *this;
    fwrite((void*)&mNx,sizeof(double),1,output);
    fwrite((void*)&mNy,sizeof(double),1,output);
    fwrite((void*)&mNz,sizeof(double),1,output);
    fwrite((void*)&mN,sizeof(double),1,output);
    
    
  for (int z = 0; z < mNz; z++) {
    for (int y = 0; y < mNy; y++) {
      for (int x = 0; x < mNx; x++) {
          nodeSP * nd   = L(x, y, z);
          nd->save(output);
      }
    }
  }
}

void latticeSP::load(FILE * input)
{
    latticeSP & L = *this;
    if(mNodes!=0) delete[] mNodes;
    fread((void*)&mNx,sizeof(double),1,input);
    fread((void*)&mNy,sizeof(double),1,input);
    fread((void*)&mNz,sizeof(double),1,input);
    fread((void*)&mN,sizeof(double),1,input);
    
    mNodes = new nodeSP[mNx*mNy*mNz];

  for (int z = 0; z < mNz; z++) {
    for (int y = 0; y < mNy; y++) {
      for (int x = 0; x < mNx; x++) {
          nodeSP * nd   = L(x, y, z);
          nd->load(input);
      }
    }
  }
}

void latticeSP::save_strings_vtk(const char * filename) {
  std::vector<double> colors;
  colors.resize(mN);

  char buff[128];
  int  x, y, z;
  int  n, k;
  int  current;
  int  c1, c2, flag1, flag2;
  int  npoints;

  snprintf(buff, 128, "%s.vtk", filename);

  FILE * output = fopen(buff, "wb");

  fprintf(output, "# vtk DataFile Version 3.0\n");
  fprintf(output, "String configuration\n");
  fprintf(output, "ASCII\n");
  fprintf(output, "\n");
  fprintf(output, "DATASET UNSTRUCTURED_GRID\n");

  latticeSP & L = *this;

  for (z = 0; z < mNz; z++) {
    for (y = 0; y < mNy; y++) {
      for (x = 0; x < mNx; x++) {
        nodeSP * nd   = L(x, y, z);
        nodeSP * ndx  = L(x + 1, y, z);
        nodeSP * ndxy = L(x + 1, y + 1, z);
        nodeSP * ndy  = L(x, y + 1, z);
        nodeSP * ndyz = L(x, y + 1, z + 1);
        nodeSP * ndz  = L(x, y, z + 1);
        nodeSP * ndxz = L(x + 1, y, z + 1);

        // Measure the XY loop
        for (k = 0; k < mN; k++) {
          colors[k] = 0.0;
        }

        measure_link(nd, ndx, colors);
        measure_link(ndx, ndxy, colors);
        measure_link(ndxy, ndy, colors);
        measure_link(ndy, nd, colors);

        flag1 = 0;

        for (k = 0; k < mN; k++) {
          if (colors[k] > 0.5 * M_PI) {
            flag1 = 1;
            c1    = k;

            break;
          }
        }

        flag2 = 0;

        for (k = 0; k < mN; k++) {
          if (colors[k] < -0.5 * M_PI) {
            flag2 = 1;
            c2    = k;

            break;
          }
        }

        if ((flag1 != 0) && (flag2 != 0)) {
          line_list.push_back(line(x + 0.5, y + 0.5, z - 0.5, x + 0.5, y + 0.5, z + 0.5, c1, c2));
        }

        // Measure the YZ loop
        for (k = 0; k < mN; k++) {
          colors[k] = 0.0;
        }

        measure_link(nd, ndy, colors);
        measure_link(ndy, ndyz, colors);
        measure_link(ndyz, ndz, colors);
        measure_link(ndz, nd, colors);

        flag1 = 0;

        for (k = 0; k < mN; k++) {
          if (colors[k] > 0.5 * M_PI) {
            flag1 = 1;
            c1    = k;

            break;
          }
        }

        flag2 = 0;

        for (k = 0; k < mN; k++) {
          if (colors[k] < -0.5 * M_PI) {
            flag2 = 1;
            c2    = k;

            break;
          }
        }

        if ((flag1 != 0) && (flag2 != 0)) {
          line_list.push_back(line(x - 0.5, y + 0.5, z + 0.5, x + 0.5, y + 0.5, z + 0.5, c1, c2));
        }

        // Measure the ZX loop
        for (k = 0; k < mN; k++) {
          colors[k] = 0.0;
        }

        measure_link(nd, ndz, colors);
        measure_link(ndz, ndxz, colors);
        measure_link(ndxz, ndx, colors);
        measure_link(ndx, nd, colors);

        flag1 = 0;

        for (k = 0; k < mN; k++) {
          if (colors[k] > 0.5 * M_PI) {
            flag1 = 1;
            c1    = k;

            break;
          }
        }

        flag2 = 0;

        for (k = 0; k < mN; k++) {
          if (colors[k] < -0.5 * M_PI) {
            flag2 = 1;
            c2    = k;

            break;
          }
        }

        if ((flag1 != 0) && (flag2 != 0)) {
          line_list.push_back(line(x + 0.5, y - 0.5, z + 0.5, x + 0.5, y + 0.5, z + 0.5, c1, c2));
        }
      }
    }
  }

  if (line_list.size() == 0) {
    return;
  }

  // Dump points
  npoints = 4 * line_list.size();

  fprintf(output, "\nPOINTS %d float\n", npoints);

  for (std::list<line>::iterator pline = line_list.begin(); pline != line_list.end(); pline++) {
    fprintf(output, "%f %f %f\n", pline -> xs, pline -> ys, pline -> zs);
    fprintf(output, "%f %f %f\n", pline -> xe, pline -> ye, pline -> ze);
    fprintf(output, "%f %f %f\n", pline -> xs + DELTA, pline -> ys + DELTA, pline -> zs + DELTA);
    fprintf(output, "%f %f %f\n", pline -> xe + DELTA, pline -> ye + DELTA, pline -> ze + DELTA);
  }

  fprintf(output, "\nCELLS %d %d\n", npoints / 2, npoints / 2 * 3);

  for (k = 0; k < npoints; k += 4) {
    fprintf(output, "2 %d %d\n", k, k + 1);
    fprintf(output, "2 %d %d\n", k + 3, k + 2);
  }

  fprintf(output, "\nCELL_TYPES %d\n", npoints / 2);

  for (k = 0; k < npoints; k += 2) {
    fprintf(output, "3\n");
  }

  fprintf(output, "\nPOINT_DATA %d\n", npoints);
  fprintf(output, "SCALARS color float\n");
  fprintf(output, "LOOKUP_TABLE default\n");

  for (std::list<line>::iterator pline = line_list.begin(); pline != line_list.end(); pline++) {
    fprintf(output, "%f\n", (float) pline -> color1);
    fprintf(output, "%f\n", (float) pline -> color1);
    fprintf(output, "%f\n", (float) pline -> color2);
    fprintf(output, "%f\n", (float) pline -> color2);
  }

  line_list.clear();
  fclose(output);
}

Box::Box(int xmin,
         int xmax,
         int ymin,
         int ymax,
         int zmin,
         int zmax) {
  mXmin = xmin;
  mYmin = ymin;
  mZmin = zmin;
  mXmax = xmax;
  mYmax = ymax;
  mZmax = zmax;
}

int Box::xmin() {
  return mXmin;
}

int Box::ymin() {
  return mYmin;
}

int Box::zmin() {
  return mZmin;
}

int Box::xmax() {
  return mXmax;
}

int Box::ymax() {
  return mYmax;
}

int Box::zmax() {
  return mZmax;
}


//~ Formatted by Jindent --- http://www.jindent.com
